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I. Introduction 


The use of uniformly or non-uniformly distributed 
pseudorandom numbers in systems simulation, statistical 
sampling experiments and analytical Monte Carlo work is by 
now well established. Numerous algorithns exist for 
producing such numbers from various distributions; for 
summaries of common techniques, see Knuth [5], Gaver and 
Thompson (2] or Ahrens and Dieter [1]. 


The user of pseudorandom numbers is usually not 
concerned with the details of the algorithm employed but 
rather with the results; a good algorithm, then, is- one 
which is fast, uses minimum computer memory and produces 
Numbers with satisfactory statistical properties. The 
search for statistically competent algorithaos for 
pseudorandom numbers has resulted in the specification of 
many so-called "exact" generators, that is those whose 
deviation from the true distribution concerned is the result 
of computer rounding errors rather than any defect in the 
method itself. Such methods for nonuniform random numbers 
are often based on the assumption that "good" uniform 
numbers are available fron an independent generator. 


Exact generators for nonuniform pseudorandom numbers 
are often quite complex and so assembly-level coding is 
often resorted to when implementing them in order to meet 
the computer time and amemory constraints on a good 
algoritha. An example is the LLRANDOM package developed at 
the Naval Postgraduate School by G.P. Learmonth and P.A.W. 
Lewis and described in [7]; it produces pseudorandom numbers 


from unifctm, normal and exponential distributions. This 
report describes an extension to the LLRANDOM package for 
Cauchy and gamma distributed numbers. 


The Cauchy gistribution has density function 


(‘) 407 <4 Jie ek ates 


and distribution function 


F(x) = 1 + 1 tan~! x. 
2 7 


While the shape of the Cauchy density resembles the nornal 
density, the tails are much heavier; in fact, Caucny <:2:.- 

- varjables have no expectation and an infinite variance. The 
density has mode at zero and often in applications the 
variates are often shifted by a location parameter T or 
scaled by nultiplying by a scale parameter S. Because of 
the heavy tails, Cauchy variates might. find applicatioa ag a 
"pathological" case in a systeus simulation study as well as 
in statistical sampling experiments for robust estimation 
techniques. See Chapter 16 of Johnson and Kotz [4] for 
further details on the Cauchy distribution. 


The gamma gistributjion with shape parameter A and scale 
parameter s has the density function ee 


A A-1 ~SX 


2) - f(x e 
(2) (x) some a.) nee 


where [(A) is Euler's gamma function 
A~1 =F 
(3) P(A) = F e dx . 


Note that [(n) = (n-1)! when n is a non-negative integer. 
If the random variable X has density (2) then 


BRXJ= Ass , 


nibh tat ea et tlt Nn Sia ie et aaa Spt 


V[XJ) * A /S/s? . 


When A = 1, X has the exponential distribution while X, 


' guitably scaled, has an asymptotically normal distribution 


as A =-> @. 


We note that if X has a 9 (A,1) distribution then X/s 
has a T(Aa,8) distribution, so we may set Ss = 1 in (2) as 
far as the generating algoritha is concerned. Tha output 
from the generator may then be appropriately scxled. 


Gamma random variables are used in a wide variety of 
applications: for analytical modeling, in reliability 
theory and for statistical testing (the chi-squared random 
variable with n degrees of freedom has the M(B )) 


- distribution). See [6] or Chapter 17 of [4]. for sore 


details. 


II. Use of the Subroutines 


This extension to LLRANDOM is composed of two 
independent IBM System/360 Assembler-coded subroutines: 
CAUCHY for Cauchy-distributed variates and GAMA for ganna 
variates. The name GAMA was chosen so,as not to conflict 
with the IBM mathematical library subprogram GAMMA which 
computes the gamma function (3). 


The basic conventions for using GAMA and CAUCHY are the 
same as in the LLRANDOM package: the invoking statements 


CALL CAUCHY -( IX, X, N) 
and CALL GAMA ( A, IX, X, ¥) 


will result in a vector X(1), «-. ,» X(N) Of Cauchy or 
rT (A,1.0) pseudorandom variates, respectively. The arg- 
ument IX is, in both cases, an integer seed to be used in 
the multiplicative congruential uniform generator employed 
by LLRANDOM. IX should be initialized just once ia the 
calling program to some positive integer value and_ should 


not be altered thereafter. 


The subroutine GAMA requires a source for normal and 
exponential deviates; these are obtained directly from the 
LLRANDOM package and so the statement "CALL OVFLOW" gust 
appear once in the calling program to initialize LLRANDOM. 
As mentioned previously, the output from GAMA must be scaled 
if the scale parameter is other than one; the following set 
of statements will thus be required to generate a vector of 
100 chi-squared variates with seven degrees of freedon: 


av 


DIMENSION X(100) 
CALL OVFLOW 
IX = 13726 


CALL GAMA ( 3.5, IX, X, 100 ) 


DO 50 I = 1,100 
X(I) = 2.0 * X(TI) 
50 CONTINUE 


END 


_ Cauchy variates are also often modified by location and 
scale parameters; since no expectations exist, however, we 
cannot refer to these parameters in terms of mean or 
variance. Subroutine CAUCHY is completely independent of 
LLRANDOM or any other subroutines so that the "CALL OVFLOW" 
statement is not necessary in this case. To use CAUCHY to 
produce a single variate C with location parameter T and 
scale parameter S we may use the statenents 


IX = 217663541 
CALL CAUCHY ( IX, C, 1) 
c=sS* C+T 


END 


Just as in LLRANDOM, linkage overhead between the 
calling program and GAMA or CAUCHY will be minimized if a 
vector of several variates is obtained at the same tine 
instead of just a single one. The gain in this case can be 
as much as 50 microseconds per variate in average generation 
time, an improvement of up to 50%. In GAMA, several 
constants must be calculated for each different value of the 
shape parameter A; these constants are saved between calls 
so that they need not be recouputed. It will thus be sore 


efficient to get several gamma variates with the same shape 
parameter before changing the A value, especially when A > 
3.0 when the setup computations are extensive (see lines 


174-246 of the program listing). 


Note that the technigues used in GAMA and CAUCHY anake 
use of so-called rejection methods so that the number of 
uniform (or exponential or normal) deviates needed to 
generate a single output deviate is random. When normal or 
exponential deviates are required by GAMA from LLRANDOM a 
vector of 10 deviates is called for; since not all of these 
may be used at the time they are generated, the balance are 
saved for the next call to GAMA. Thus, reinitializing the 
seed IX to its original value will not in general result in 
an exact repetition of the generated gamma sequence since 
the first few deviates will use the old normal or 
exponential deviates from the previous sequence. To achieve 
an exact repetition, the qenerator must be forced to repeat 
the initialization computations for the desired A value; at 
this time any remaining variates from LLRANDOM = are 
discarded. An example of this might be 


DIMENSION G (100) 
CALL OVFLOW 
IX = 12345 


CALL GAMA ( A, IX, G, 100 ) 


Cc REINITIALIZE GAMMA SEQUENCE 
CALL GAMA ( 1.0, IX, G, 1 ) 
IX = 12345 


CALL GAMA ( A, IX, Gy, 100 ) 


CAUCHY requires 552 bytes and, as mentioned previously, 
is completely independent of any other subprograms. CAUCHY 
uses the LLRANDOM aultiplicative congruential uniforn 


generator but this is coded in line when needed so as to 
preserve CAUCHY's independence. The average generation tine 
per variate for subroutine CAUCHY on a System/360 Model 67 
under OS/MVI was 67.5 microseconds when variates were 
generated in vectors of 100. The generation of variates one 
at a time increased the average time to 119.3 microseconds 


per variate. 


Subroutine GAMA itself uses only 1988 bytes of memory 
but since it calls on LLRANDOM the total core reguirement is 


9342 bytes: 


GAMA 1988 bytes 
LLRANDOM 6189 bytes 
Reguired IBM Functions 1165 bytes 
Total 9342 bytes 


Timing the gamma generator on a System/360 Model 67 was 
carried out using the TIME macro; Table 1 summarizes the 
observed times as a function of the shape parameter, A. 
Note that since special methods are employed when A is 0.5, 
1.0, 1.5, 2.0 or 3.0, the times in these cases are 
considerably shorter than tiaes for nearby values of A. 


Shape Parameter Algorithm © Vector of 100 = Single 


| A | - ' "* Vardates -  - ‘ Variate 

0.1 Gs 324.0 364.0 ; 

0.3 | 6S. Sb TO. 402.5 

0.5 GA 70 8099 

0.8 GS 439.8 551.2 

0.9 Gs 459.0 > 649.6 

1.0 GA | 68.7 158.9 

1.2 GF o  . 300.1 385.0 

1.4 SR SDE eco 

1.5 GA é 1 Te 215.8 

1.8 GF 343.6 = 3908 

2.0 Bk. ge: 142.5 203.6 

261 GF 396.1 . 450.8 

2.5 GF 434.7 468.5 

2.9 GF 4445... 496.6 

3.0 GA: 206.7 eee cis 
3.1 GO 8=tt—<“<iéiS NS 8 ; 

3.5 Go 3, S362 «373.4 
4.0 Go 332.4 420.7 ; 

5.0 Go S07.7 8682 

8.0 GO 293.1 371.3 

10.0 ele) 289.4 ; 312.5 

20.0 GO 238.2  - 321.6 

50.0 GO : 197.7 : 284.2 

100.0 - “66 178.4 . 220.0 

1000.0 “BOs © 166.7. 177.0 

10000.0 | fe) 3 136.4 169.8 

100000.0 GO 482.5 235.8 


Table 1. Average generation times (microseconds) for gamma 
variates using subroutine GAMA. 


III. Description of the Algorithas 


This section describes the actual algorithms used in 
CAUCHY and GAMA. An understanding of the algorithms is not 
necessary for use of the package but they are set forth here 
both in the interest of completeness and in an effort to 
document’ the programs more fully. A single algoritha 
suffices for the Cauchy generator while GAMA uses one of 
four algorithms, depending on the value of A. 


In the descriptions which follow, the letters U, N and 
E (with or without affixes) represent uniform, standard 
normal and unit exponential pseudorandon deviates, 
‘respectively. The phrase "Generate U" implies that U is the 
next sequential uniform variate in the linear congruential 
sequence; these variates are generated as needed by using 
the same aultiplicative congruential scheme as used in 
LLRANDOM. The phrases "Generate N" or "Generate E" imply 
that normal or exponential variates are to be obtained by 
linking directly to LLRANDOM. 


A. Cauchy Generator 


The Cauchy generator is a combination decomposition- 
rejection method (see Knuth [5)). The Cauchy density is 
decomposed, as in Figure 1, into three subdensities: a 
uniform density between 0 and 1 (fi). a wedge-shaped density 


eT a long tailed density (f)- 


The uniform density ‘, is sampled with probability 


1/73; in this case a uniform(0,1) variate is returned. The 


density 3 is dealt with by using MarsSaglia's almost-linear 


Ae Fhe aac 


density algoritha, just as in Knuth's Algoritha L (5). The 
density f, is sampled with probability 1/2 - 1/1 . The tail 
density ‘. is sampled by a rejection method with probability 


1/2. The majorcizing density for - is g(x) = 1 / x2, which 


is the density of the reciprocal of a uniform (0,1) variate. 


Algoritha cC below uses the fact that in the priae 
modulus congruential random nurber generator used in 
LLRANDOS the low order bits are uniforaly distributed so 
that b. and :. select the proper sub-distribution in Step 1. 


This will not in general be the case for other congruential 


pseudo-random nuaber generators. 


f(x) 


Figure 1. Decoaposition of the Cauchy Density Function. 


en Re ee ee ee 


ee ee otis wiatsee setae 


Algorithm ¢C. Cauchy variates. 


1. (Select subdensity) Generate U, setting aside the two 
low order bits a and If b. = 1, go to Step 6. 


2s (Sample box) If 0 S$ 0.6366197724 = 2/7, generate a 


* * 
new variate U , set x = U and go to Step 8. 


3. (Sample wedge) Generate new variates We and ee If Ue 


> U_, exchange UV andU. Setx=JU. 
2 1 2 1 
4. (Easy vrejection) If a S$ 0,8284271247 = 2/2 - 2, go 
to Step 6. 
5. (Hard rejection Tf vu -U $s 1 = x? 2/2 - 2 o to 
( 3 ) "2 1 qe ( 2/2 ye J 
Step 8, otherwise go back to Step 3. 


$6. (Sample tail) Set x= 1/ U. 
* ea 
Ve (Tail rejection) Generate a new variate 0. If JU s 


on go to Step 8, otherwise generate a new U and go 
back to Step 6. 


8. (Random sign) pa a, = 1 set x = - x. Deliver x as the 


generated deviate. 


It should be noted that there are several other methods 
for generating Cauchy variates: the ratio of independent 
standard normal deviates has the Cauchy distribution, as 


does the quantity 


X=t 1? ee ecg P) 
an{™ ( 1 


where U is uniform (0,1). These methods are both substant- 
dally slower than algorithm C, but another new method has an 


_ average time comparable to Algorithm C and is much easier to> 
prograa, This second method requires an average of 2,55 
uniform random variates per Cauchy variate (as compared with 
2.47 for algorithm C) and it needs about 69 microseconds per 
variate on the Systen/360 Model 67. It is possible, 
however, that Algorithm CR will be better than algoritha ¢ 
‘in soue other imapleaentation, 


; The method is essentially the technique devised by von 
Neumann to generate a random variate sin J, where U is 
uniform between 0 and 2n. Such variates are used in the 
polar method for generating normal random variables (8]. It 
| does not seen to have been recognized that the method also 
generates tan U, which is the required Cauchy variate. 


Algoritha CR. Cauchy variates, ratio method, 
1. (Get uniforas) Generate 2 and eas Set - = 2 uf - 1 
and Y #20 <1. 
2 2 
2. (Rejection test) If be + pr > 1 go back to Step 1. 


-3. (Take ratio) . Deliver x = Y 7%) 


Be. Gampa Generator GS: A < 1.0 


This method is due to Ahrens and is set forth in [1]. 
It is applicable only to values of A less than one and is 
markedly superior in execution time to the method of Johnk 
{3], which is the usual technique for generating variates of 
this type. | ; . 


The method is a rejection method employing two 
different tests, one of which is chosen at random for any 
given variate; the power transform of a uniform(0,1) 


s/h ; 
variate, U m . is tested in the region 0 < xX 1, while a 


Suitable exponential, E, is tested when x > 1. The 
advantage of this method lies in the limited use of the 
library subprograms for the exponential and logarithn; 
average times range from 300 to 400 microseconds as compared 
with 600 to 800 for Johnk's method. Purther discussion and 
proofs may be found in {1}. 


Algorithm G§. Gamma variates, A < 1.0. 
1. (Select rejection test) Generate U and generate £ and 
set P= e+A UU. (Note that "e® is the base of the 
CC) 


natural logarithas.) If P $ 1 go to Step 2, otherwise 


{ 
go to Step 3. 


fa 
2. (Small x_ test) Set x = P e « Ifx S$ E, deliver x, 


otherwise go back to Step 1. 
3. (Large x _ test) Set x =- ln [ } {e+ An~ P} } IF 
e 


(1 - A) In x $ E, deliver x, otherwise go back to Step 
1. 


A thus-—far unpublished method devised by Professor G.S. 
Fishman of North Carolina University was communicated to the 
authors in private correspondence. It is valid for any A > 
1.0 but its efficiency in terms of average time goes down as 


Vk so it is applied in GAMA only in the range where it is 


superior to the Dieter-Ahrens method GO described below. 


The method is a_ rejection method based on the following 
_ theorea. 


Theorem Let G be a uniform (0,1) random variable and let £ . 
be an exponential random variable with mean A. Let 


A= <x (1-474) = (A=1) 
g(x) = ¢ xq. e ‘ 


If g(E) 2 U, then E has conditionally the gamma distribution 
with shape parameter A, i.e. 


A-1 -x 
; zt x | US$ g(B) ) =x : er . 


Proof: - ; 


P, . = x A 
Unconditionally, E has density h(x). = } e id ° 


Therefore, 


(4) (x | USg Bix) Fe fpsa tat peash . 
e 


Now since U is uniformly distributed, , 


Pr (USg (E) JE=x} = g(x) 


as long as 0 < g(x) < 1; that this is true for every x > 0 
may be readily verified by elementary calculus. Therefore, 


(5) «Pe (USg(Z)) = EL Pr(USg(E)1E} ] 
= I, g(x) h(x) az" 
ye A ee 
= C{A) 


Thus, in view of (4), ; é 


a 

H. 
le. 
a 
a 
ig 


£ , XI 0S¢ (2) = Bix) g(sh 


A-1 -x 


Dear ac ¢ 4 


The efficiency of the generator is governed by the 
probability that a given variate will pass the rejection 
test, U Ss g(E)3 from (5) it will be seen that this 
probability is just C(A). When A is large we have fron 
Stirling's approximation that C(A) = rose , so that the 

a 
method becomes more inefficient with increasing A, as noted 


é 


above. 


A slight modification to the method suggested by the 
theores improves the efficiency slightly and we obtain 


Algorithm GF. Gamma variates, 1.0 < A < 3.0. 


Se (Generate exponentials) Generate two independent 


exponential variates, , and zs 
2. (Rejection test) If B < (A-1) (Z, - la e. - 1) then 
go back to Step 1. 


3. (Acceptance) Deliver x = A oF 


D. Gamma Generator G05 A 2 3:9 


This method was originally developed by Dieter and 
Aurens and is fully described in [1] together with several 
other gamma generation techniques. Algorithm GO does not 


suffer the usual drawback of growing less efficient in 
generation tine with increasing A; in fact, the amethod is 
more efficient for larger A values. : 


The basic idea here is to take advantage of “the 
asymptotic nocmality of the gamma distribution by doing nost 
of the sampling from a norawal distribution; the right hand 
tail is sampled, when necessary, using a rejection method 
with the exponential distribution. The method can be 
applied to values of A greater than 2.533, but it is not as 
efficient as Fishman's technique for A < 3.0. 


As mentioned previously, this algorithm requires the 
computation of several constants which depend only on A ana 
which may be saved between calls; these calculations are 
déscribed in step 0 of the specification below. Further 
discussion, illustrations and proofs are given in [1]; the 

“version of GO here differs in a few minor details from the 
original Dieter and Ahrens technique. 


0. (Calculate constants) Coapute: 
mea- 1; 


s2 = + ie s = Ss?; 
x de 


ad = /ész; “pede ms. 


w=s2/an- i; v = 282 / (m YA )5 


c= b¢il1n s_4 - 2m - 3.7203285. 


1. (Select normal/exponential) Generate JU. If us 
0.0095722652 go to Step 7. 

2. (Normal sampling) Generate N and set x = SN + a. 

3. (Check trial value) If x < 0 or x > b go back to Step 2, 
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otherwise generate a new variate U and set S = N2@ / 2. 
If N > 0 go to Step 5. 

4, (Left-hand rejection) If 0 < 14 S (vN - w) go to Step 
9, otherwise go to Step 6. 

5. (Right-hand rejection) Ifu < 1-+wS go to Step 9. 

6. (Final normal rejection) If ln U < mw ln x +ma-x+ § 


go to Step 9; otherwise go back to step 1. 
7. (Exponential) Generate a. and os and set x = b(1¢8, 7/4) s 
8. (Exponential rejection) If » (x - ln x ) + c> E. go 
back to Step 1. 


9. (End) Deliver x as the gamma variate. 


E. Ad Hoc Gamma Generators 


This set of algorithuas is based on the well-known fact 
that the sum of independent gamma variates with shape 
parameters a; and ‘ and equal scale parameters has the 
gamma distribution with shape parameter ig A, and scale 
parameter equal to that of the summands. We may thus 


generate a gamma variate with integer shape parameter K by 
taking the sum of K independent exponentials. This will be 
more efficient than the previously discussed methods 
(Algorithwus GF and GO) for moderate values of K; for the 
System/360 we take K < 3 to apply this ad hoc technique. 


An obvious extension to this method is to allow for 
half-integral values of A by making use of the fact that the 
square of a standard normal random variable has the 
chi-squared distribution with one degree of freedom, i.e. 
N2/2 has the ganna distribution with unit scale parameter 
and A # 0.5. We use this extension for A = 0.5 or 1.5. 


The resulting algorithm is then 


Algoriths GA. Gamma variates, integral or half<integral 


AD 


2. 


4, 


botacte dt 6 babe so uaks 


shape parameter A. 


' (Pind K) Set K = (A), where (A} denotes the integral 


part of Aw Set X= 0. Tf A-K= 0.5 set L = 1; if A- 
K = 0.0 set L = 03; otherwise Stop. (If the algoritha : 
stops, an incorrect A value has been used.) 

(Generate exponentials) If K = 0 go to Step 3, 
otherwise generate K exponentials Ee cee hare and set 


X & E ¢t eee + B e 
1 K 


(Generate normal) If L = 0 go to Step & otherwise 
generate N and set X = X + N2/2. 
(Deliver X) X is the desired variate. 


IV. Summary and Conaents 


. This work provides a convenient and useful extension to 
the LLRANDOM package, especially for users interested in 
. statistical. and reliability theory applications of digital 
Simulation. The combination of the most efficient known 
gamma generation techniques with the new Cauchy method gives 

exceptionally good time characteristics at some cost in 


computer memory utilization. 


The work may be extended at once to the generation of 
several other types of random variables. For example, the 
beta distribution with parameters A and B may he sampled by 
taking gamma variates x and x. with respective shape 


parameters A and B and delivering. 


Zz X Xx 
* ~ 1 : 2 


as a beta variate. In this case considerable overhead in 
GAMA can result froa shifting the shape parameter back and 
forth between A and B; for this reason obtaining vectors of 
gamma variates X and x. is recommended, as in the following 


exaaple: 


DIMENSION X1(50), X2(50), 2(50) 


CALL GAMA ( A, IX, X1, 50 ) 

CALL GAMA ( B, IX, X2, 50 ) 

DO 405 I = 1,50 

Z(I) = X1(T) Z ( XU(I) + X2(T) ) 
405 CONTINUE 


END 


The t-Distrpibution may be sampled as the ratio of a 
standard normal and an independent chi-squared random 
variate, while the F-Distripution may be obtained by tcking 
the ratio of two independent chi-squared variates divided by 
their respective degrees of freedom. (See pages 4 and 5 for 
an example of the generation of chi-squared variates.) 


EE ea ee ee ee nr ee eee 
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